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SUMMARY 


Dynamic systems that were once controlled by analog circuits are now controlled by digital 
computers. Presented is a comparison of the digital controllers presently used with magnetic 
suspension and balance systems. The overall responses of the systems are compared using a 
computer simulation of the magnetic suspension and balance system and the digital controllers. The 
comparisons include responses to both simulated force and position inputs. A preferred digital 
controller is determined from the simulated responses. 
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1. INTRODUCTION 


The first recorded use of an actively stabilized magnetic suspension system was at the 
University of Virginia, USA, in 1937 (Ref. 1). Such systems are now finding many uses, 
including the suspension of models in wind tunnels. 

Magnetic suspension of a model in a wind tunnel was first achieved in 1957 by researchers at 
the Office National d’Etudes et de Recherchcs Aerospatialcs (ONERA), France (Ref. 2). The 
ONERA system controlled models in five dcgrces-of-frcedoin in test sections up to 30 cm in 
diameter. So far as is known, 17 wind tunnel magnetic suspension systems have been built since 
then, with six now in operation (Ref. 3, 4). 

All wind tunnel magnetic suspension and balance systems (MSBSs) use controlled dc 
electromagnets acting on a suspended body containing a ferromagnetic material. With this 
approach, stabilization of the position and attitude of the suspended body requires feedback 
controllers. Early control systems used analog circuits, each individually designed for a 
particular system. Performance was restricted by practical limits on complexity and adjustment 
of the controller parameters, and stability of the analog elements. With the development of 
digital computers, digital control became possible, promising many advantages. 

One advantage of a digital controller is that it requires less hardware than an analog 
controller. A digital controller uses digital-to-analog (DAC) and analog- to-digital converters 
(ADC) for communication between the computer and the MSBS. The control strategy is 
written in software and is easily modified to improve control techniques, either as better 
computer systems become available, or the MSBS changes. With a digital controller 


1 


the possibilities of controllers are limitless and the great flexibility of software allows complex 
control strategies. 

Of the six known wind tunnels using magnetic suspension and balance systems, two are at 
NASA Langley Research Center in the USA. The others arc at Oxford University and the 
University of Southampton in England, the National Aerospace Laboratory (NAL) in Japan, 
and The Central Aero- Hydrodynamics Institute (TsAGI) in the Soviet Union. 

All of the existing MSBSs arc fitted to relatively small wind tunnels. The largest system, 
which is in the Soviet Union, installed in a 40 x 60 cm test section and is used for low speed 
aerodynamic testing (Ref. 5). Both of the MSBSs in the USA are fitted to low speed 
atmospheric fan-driven open-return tunnels. One of the USA MSBS wind tunnels has a 15 cm 
diameter octagonal test sect ion. The other, known as the Langley 13 inch MSBS, has a 26.7 x 
31.8 cm octagonal test section and is used on a fairly regular basis for low speed aerodynamic 
testing. The MSBS at Oxford is (1 f ted to a 12 x 12 cm hypersonic tunnel. The most highly 
developed MSBS is at the University of Southampton. The Southampton system is fitted to an 
18 cm octagonal test section and is used for dynamic as well as static aerodynamic testing. The 
newest MSBS is the NAL system which is fitted to the 10 x 10 cm transonic test section of their 
Pilot Cryogenic Wind Tunnel. 

Of these six MSBSs, only three arc digitally controlled. These are the University of 
Southampton MSBS, the NASA Langley 13 inch MSBS, and the NAL MSBS. The 
Southampton MSBS digitally controls 10 electromagnets using a minicomputer to maintain 
control of the model in six degrees of freedom. The NASA Langley 13 inch system has only 5 
electromagnets controlling five degrees of freedom. The NAL system controls five degrees-of- 
freedom using 10 electromagnets. 

Table I gives a complete listing of the existing MSBS wind tunnels. 
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2. GOVERNING EQUATIONS 


The principles of ail MSBS can be understood by studying a single degrcc-of- freedom 
system. Figure 1 shows a simple single degree-of-freedom MSBS consisting of a dc 
electromagnet and a suspended magnetic body. The suspended body must contain some 
ferromagnetic material. The electromagnetic field from the coil produces a magnetic force which 
attracts the suspended body to the coil. Gravity acts to pull the suspended body away from the 
coil. If the current in the coil increases, the magnetic force of attraction also increases. 



9 


Figure 1. Schematic of single degree-ofTreedom MSBS. 
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Force of attraction, 


I 

As shown in figure 2, for constant coil current, the magnetic force attracting the body 
decreases as the separation distance, x, increases. This decrease in the magnetic force attracting 
the body as the separation distance increases makes this system inherently unstable. Because 
this system is inherently unstable, a feedback control system is required to regulate the coil 
current. The control system must increase the current when the separation distance increases 
and reduce the current when the separation decreases. Stable suspension of the body is possible 
through proper regulation of the current by the controller. 



Figure 2. Magnetic force - distance characteristics at constant current. 
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2. 1 Dynamics of the Suspended Body 


The equation of motion for the suspended body is derived from Newton’s second law of 
motion. 


mx^£F 


Neglecting bouyancy, there are four forces acting on the suspended body in a single degree~of- 
freedom system as shown in figure 1. These forces are gravity, the magnetic force produced by 
the coil, a damping force, and any external force acting on the body. Taking positive x in the 
direction of gravity, the equation of motion for the body is: 


m x = F g - F^(x, i) - + f ( 2 -*) 

In equation 2.1, F g is the weight of the body, F A is the magnetic force exerted on the body 
by the coil, Fp is the damping force acting on the body, and f is an external force. 
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The magnetic force F^, is usually nonlinear. It is a function of the current in the coil and 
the position of the suspended body. Figures 2 and 3 show how this magnetic force varies with 
coil current and position of the body. The variation in force with x and i may be linearized by 
limiting the motion of the body and the current in the coil to small variations around their 
equilibrium values. (Ref. 6, 7) 

Let i(t) = i 0 +<5i(t) where i 0 is a constant current and 5i(t) is a small time-dependent 
variation in current around i 0 . Let x(t) = x 0 -f£x(t) where x 0 is an equilibrium position and 
6x(t) is a small variation in position around x 0 . Therefore: 


^ *o) + j\)l 


*0 


MO + |(f a ) 


X 0 . 'o 


6i(t) + higher order terms 


( 2 . 2 ) 


F(x 0 , i 0 ) is the magnetic force of attraction caused by the current i 0 with the body at an 
equilibrium point x 0 . The partial derivatives of F^ are the slopes of the force curves for 
constant current and constant position. Under equilibrium conditions, l (x 0 , i 0 ) is the magnetic 
force required to exactly balance the gravitational force acting on the body and any external 
forces which arc constant. Therefore: 


H x oi *o) — Lg — mg + Constant 

9 d 

For small variations in current and position, let ^(F^) x i 0 = anc ^ 5i^A^x 0 , i 0 = ^i* 

Equation 2.2 is further simplified by neglecting as insignificant the higher order terms. As 
shown in figures 2 and 3, the slo|>es of the curves are such that K x is negative and Kj is 
positive. These force constants can be considered to represent the spring-like stiffness of the 
system. These force constants can !>e determined experimentally for a given equilibrium current 
and |K>sition. 
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Equation 2.1 takes account of the damping forces acting on the body caused by both 
aerodynamic (viscous) and eddy current damping. These damping forces are assumed to be 
velocity dependent. The eddy current damping is usually very small and can be ignored. 
However, the aerodynamic damping can be large, especially for wind tunnel testing. The 
damping term has a negative sign because the damping force always opposes the motion. A 
motion in the positive direction produces a damping force in the negative direction and a motion 
in the negative direction produces a damping force in the positive direction. With small 
variations in position, the damping force becomes: 

F D = C X 

The linearized equation of motion for the sus|>ended !>ody about an equilibrium point is: 

in £x(t) = K x Sx( t) - Kj S\(t) - C «x(t) + f (2.3) 

In this equation 2.3, f is the change in external force. 

2.2 Governing Equation of the Magnetic Coil 

The governing equation of the electromagnetic coil is the sum of the voltage drop across the 
coil resistance and the voltage across the electromagnetic coil. 

V(t) = i(t) R + ^i(t) l) = i( l ) + I' (2-4) 

Where V(t) is voltage, i(t) is current, L is inductance, and R is resistance. 

In addition to being a function of the geometry of the coil, the inductance of the coil is a 
function of the suspended objects position, L = L(x). The time rate of change of the inductance 
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can I* simplified by invoking the chain rule, ^(L(x)) = £{l) &(x(t)} Substituting this in 


equation 2 A gives: 


V(t) = i(t) R + L £(i(l)) + i(t) ±(l) < f l (x(t)) (2- 

This velocity, ^( x ( t )). > s caused by changes in the inductance L, resulting from the motion of 
the body. This velocity is not related to a change in coil current. (Iter. 8) 

One method of linearizing equation 2.5 is to assume V(t), i(t), and x(t) are allowed only 
small variations around some equilibrium points as assumed in the equation of motion for the 
suspended body. For small variations, V(t) = V 0 +*V(t), i(t) = i„+*i(l), and x(t) = x 0 +5x(t). 
Substitution of these expressions into equation 2.5 gives: 


v 0 + «v(t) = (i„ + 6i(t))lt + L ^ t (i„ + *i(t)) + (i 0 + *i(t))£(L) ^(x o +*x(0) 

V 0 + «V(t) = i 0 R + «(0 R + b + '0 «(t) + «i(t) ^(L) 6x(t) (2.6) 


Since V 0 = i 0 R t this becomes: 


tfV(t) = ^i(t) R + L ^(«i(t)) + i 0 ^(L) M(l) + «(0 ^(b) ( 2 - 7 ) 

If <5i(t) and <5x(t) are very small, then their product is even smaller and can be neglected as 
insignificant. Equation 2.7 is further simplified by letting 'o^“(b)lx 0 = K c because j^(L) is a 
constant slope for small changes in position as shown in figure A. 
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Coil inductance 



Figure 4. Induction - body position characteristics at constant current. 

Therefore the linearized approximation of equation 2.4 is: 

SV( t) = 5i(t) R + L ^L(«i(t)) + K c Sx( t) (2-8) 

2.3 Single Dcgrec-of-Frecdorn MSBS Tran sfer Function 

The system differential equations for small variations arc equations 2.3 and 2.8. 


m 6x(t) = K x <$x(t) - Kj 5i(t) - C £x(t) 4- f (2*3) 

tfV(t) = Si(t) R + L + Kc ^x(t) (2-8) 


Assuming the initial conditions are zero, 


these equations transfer to the Laplacian S-domain as: 


m S 2 AX = K x AX - Kj AI - C SAX + f 


AV - AI R + L 5AI + K c SAX 


S 2 + C S - K x ^ AX = - Kj AI + f 
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AV = Al (R + L s) + K c 5AX 



AX(m S 2 + C S - K x ) = -KjAI + f 


(2.9a, b) 


AI = 


AV - K C S AX 

It + L S 


AX (s 2 + gs-fe) = 


-K. 


m 


-AI + m 


AI = 


AV 

_ IT 


-jr SAX 


( l + fe s ) 


AX 


K i K } K c 

(5 2 + ^ 5 - fe) = ( ”7 -X AV + - m f -r 

v ' ( 1 + H ( , + it 5 ) 


5AX + 1 




AV + j[j 


Combining equations 2.9a and 2.9b gives the transfer function of this single degree-of-freedom 
system (in control nomenclature, this is referred to as the plant transfer function): 



(2.10a) 
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AX = 


m L 


AV + 


R ( 
mTA 


1 + 


k s ) 


S 3 + S 2 


fR , C\ , JcR . Kx . K i M . R_K> 
\L + m / + I m L m m L J TTm 


(2.10b) 


A block diagram of this plant is shown in figure 5. 

For the system with no change in external force inputs, f=0, there are three poles. The 
poles are coupled as seen in equation 2.10a. The pole located at -R/L is the lag time created by 
the power supply and electromagnetic coil. 

The other two poles depend on the constants associated with the MSBS and the lag time. 
Typically these two poles are paired in the complex plane with a pole to the right and a pole to 
the left of the imaginary axis. The positive pole causes the system to be unstable. 



Figure 5. Block diagram of MSBS plant. 


12 








Imaginary 

A 




Real 


Figure 6. Location of MSBS plant poles in the complex plane. 

Figure 6 shows the pole locations of a linearized MSBS plant. By observing the effects the 
system parameters have on the pole locations, it is possible to modify the design of a MSBS to 
position the poles. 

The resistance of the coil, R, has a large influence on the location of the pole Pj. Increasing 
R will move Pj to the left in the complex plane. Increasing R will also move P 2 slightly to the 
left and P 3 slightly to the right. 

The inductance of the coil, L, also has a large influence on the location of V v Increasing L 
moves P, to the right in the complex plane. Increasing L will also move P 2 slightly to the right 
and P 3 slightly to the left. 

The negative value of K x is the primary reason for the instability of a MSBS. Increasing 
the negative value of K x causes the poles and P2 to move to the left while moving pole P3 
to the right. 
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Increasing the damping coefficient, C, moves the poles Pj and P 3 to the left and P 2 to the 
right. This increased aerodynamic damping usually increases the stability of the MSBS. 

Another parameter often available during the design of a MSBS is the mass, m, of the 
suspended body. Increasing the mass moves the poles P 1 and P 2 to the right, and pole P 3 to 
the left. 

The constants 1C and K c will shift the poles in the same directions. Increasing Kj or K c 
causes pole Pj^ to move left, and poles P 2 and P 3 to move right. 

2.4 State Space Representation 

The system differential equations with a small input force disturbance, f, are: 


6x(t) = ^ 6x(l) - nMi(l.)'£**( l ) + iTi 
6i(t) = l5V(t) - S 6i(t) - ^ <5x(t) 


By choosing the state variables as 6x, 6x, and <*>i , the state-space form is: 


■ 





- 


- 

— 

Sx 


0 

1 0 


bx 


0 

0 

6x 

- 

K x 

m 

-C K i 
m iri 


bk 

+ 

0 

1 

m 



0 

-k c .a 
17 T 


b\ 


1 

1, 

0 j 

_ 










6V 

r 


y - I i o o J 


6x 

6i 

Si 
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For f=0 the system equations are: 


6x(t) = 5x(t) - ^ 6 i(t) - £x(t) 5i(t) = f/V(t) - ^ 5i(t) - ^ 6x(t) 


K r 


ami the state-space representation is: 


- 


- 


— 


' 


" 

<5x 


0 

1 

0 


fix 


0 

6x 

= 

K x 

TTT 

- C 
m 

- Kj 

ttt 


6x 

+ 

0 

si 


0 

- K c 

l 

- R 
17 


6 i 


l 

l 



( 2 . 12 ) 


y = [ 1 0 0 


Sx 

Sk 

6\ 


This state-space representation can he shown to I >o controllable and observable. Because 
this system is controllable and observable, state-space control laws can be used to control the 
system. With a state-space controller the poles of the controlled system can be positioned at 
any desired location in the complex plane. 
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3. MAGNETIC SUSPENSION AND BALANCE SYSTEM CONTROLLERS 


The typical MSBS is a multiple degree-of-frecdom system using, as a minimum, one 
electromagnetic coil for each degree-of-f'eedoni controlled. 



vertical 

force 


information 


Pigun- 7. MSBS control loop. 


The controller in an MSBS used with a wind tunnel must stabilize and control the axial, 
lateral, and heave (x, y, z) positions and the roll, pitch, and yaw (<£, 0, xp) orientations of the 
suspended model (although roll is often left open-loop). This requires continuous adjustment of 
the currents in the electromagnetic coils. The adjustments of the coil currents must modify the 
attraction force curve in figure 3 to that shown in figure 8 below. (Ref. 9) 


Hi 



Force of attraction, 


< 

U. 


Separation distance, x 

Figure 8. Magnetic force - distance characteristics as modified by the controller. 

From the plant transfer function given in equation 2.10, it can be seen that the system is 
inherently unstable. A position feedback is insufficient to achieve stability, therefore some form 
of rate information is necessary (Ref. 10). 

Because position information is usually available, the traditional approach for an MSBS 
controller is to generate limited rate information ([Kxsition derivatives) using analog phase- 
advance controllers, proportional-derivative controllers, or a proportional-integral-derivative 
controllers, often combined with error integrators to minimize steady-state errors. The 
controller is located either in the forward path or the feedback path. 

3.1 Phase Advance Controller 

The standard form of a phase-advance controller is: 

input -► A ^ output 
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(3.1) 






Where A and B are the phase-advance time constants and the ratio of A/B is the high-low 
frequency gain. 

A single phase-advance can be adequate for some systems, although two or more are usually 
combined in series. The values of A and B would depend on the pole locations of an MSBS 
plant and the desired system performance. 

A single phase-advance has one pole and one zero. The pole and zero of the phase-advance 
controller should be located so they affect the stability of the MSBS plant. The idea is to choose 
a zero for the phase-advance which will make the system stable. Figure 9 shows the 
modifications that a phase-advance makes to the root locus, giving the system a stability range. 
The actual location of the pole and zero will be based on the plant poles and the desired system 


performance. 



Real 


Figure 9. Root locus of MSBS with phase-advance controller. 


i ! 
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3.2 Proportional-Integral-Derivative Controller 


The standard form of a proportional-integral-derivative (P I D) controller is: 


input 



+ K d J + 



-► output 


(3.2) 


This controller will have a pole located at the origin of the complex plane and two zeros to 
the left of the imaginary axis. Again, the location of the zeros can be selected to provide a 
range of stability for the system. Figure 10 shows how a P I D controller modifies the root locus 
of the MSBS plant. 



Figure 10. Root locus of MSBS with P I D controller. 
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4. CONTROLLERS 


4.1 Development of Digital Controllers for Wind Tunnel MSRSs 

The use of digital controllers in MSBSs allows an infinite number of possibilities for 
controllers. The first attempts at using digital controls were simply to simulate existing analog 
control systems. The approach of digitally simulating the analog controller can be simple or 
complex as shown for the following MSBS digital control systems. The sections which follow 
present a chronological history of the development of digital controllers for wind tunnels. 

4.1.1 Oxford, England: 

The development of digital control systems for an MSBS started in 1971 at Oxford 
University. The Oxford MSBS controller was implemented with conventional circuitry, using 
analog sample-and-hold stages. Discrete- time control was necessary due to the use of a scanning 
TV system for position detection of a small sphere (Ref. 11). Three degrees^of-freedom, the 
horizontal and vertical position, were controlled in the MSBS. Although it did not use a true 
digital controller, the work is noteworthy since it was founded on the same theoretical basis as 
later digital controllers. Eurthermorc, the system required a formidably complex piece of 
circuitry. 

The control algorithm is derived from a ^-transformation of a phase-ad vance controller. The 
phase-eld vance transfer function expands in the r-dornain as: 

1 “h a_ ^ z -}“ a_ 2 z 
1 ~h l). 1 z -f 1)_2 <2 ^ 
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rewritten as a difference equation, this transfer function is: 


V k = K [ f k + a -l C k-. + a -2 <k-2 " b -> V k-. - b -2 V k- 2 ] 

The controller was located in the forward loop of the system. The output, V^, is based on 
the previous and twice previous command signals and the present, previous, and twice previous 
error signals. The system used 100 control cycles per second. 

This system was later developed to include an integrator in the forward path with 
combinations of phase-advance controllers (Ref. 12). 

4.1.2 MIT, United States: 

The next developments occurred at MIT in 1976, when the theoretical application of full 
digital controls to the MSBS was studied (Ref. 13). MIT developed a hybrid simulation of an 
MSBS using a microcomputer and an analog computer. A one degree-of- freedom demonstration 
system was digitally controlled using a z - transformation of a triple phase-advance controller on 
an INTEL 8080 microprocessor. The single degree-of- freedom triple phase-advance controller 
had the following form: 

(z-a 3 )(r-ci 2 )(> ai ) 

)_ 

The researchers at MIT gave guidelines for the computing power required for a full MSBS 
system. However, financial support could not be obtained for further development of this 
system and the work was dropped. 
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4.1.3 Southampton, England: 


In 1981, researchers at the University of Southampton developed a two degrcc-of-frccdom 
digital controller for their MSBS (Ref. 14). Initially only vertical translation and pitch rotation 
were controlled by digitally simulating an analog dual phase-advance controller using a PDP- 
11/34 computer. The Southampton system placed the controller in the feedback path of the 
circuit and an error integrator in the forward path. The dual phase-advance transfer function is: 


V _ 

t ~ 


(1+ nA S) 2 
(1+ B Sf 


The Southampton digital algorithm is derived from a difference equation approximation of 
the controller transfer function. The transfer function is split into four blocks where the third 
and fourth blocks are the same as blocks one and two. The phase-advance time constants A and 
B are equal. The n is a constant to obtain the desired high/low frequency gain, nA/B, for the 
phase-advance controller when A and B arc equal. 


c 



y -♦ 


[ 1+nA 5] 2 -V'-c' 



[ 1+ nA S] 4 -* V 


The first two blocks were originally approximated as follows: 


A 


Ay _ 
T " 


£ k " y k-i 



= y k + ,,A 


Ay 

~T 


where Ay = y k - y kl 

giving y k = (J) f k + ?!;-! V *k = ) y k " ( V) y k-i 
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fjl rp T'f’D A A 

^ ai= A’ a2Z= “Jv~’ a 3“ — T*"’ ant * a 4 = ^ then these equations can be resolved into 
difference equations where: 


y k “ a i f k + a 2 y k-i 


V 'k = a 3>k + a ^k-. 


and 


y k-i - a l f k-i + a 2 y k-2 V 'k-1 = a ayk-i + a ^k-2 

then combined: 

V k = a 2 V 'k-i + a 3 a i f k + a 4 a l f k-l (4.1a) 

Also from the third and fourth block: 

V k = a 2 V k-i + a 3 a i f ' k + a 4 a ifVi (4-lb) 

Assuming V —t 1 and combining equations 4.1a and 4.1b in series is then a difference 
approximation of a dual phase-advance. The values of the constants A, n, and T used were 
different for the two degrecs-of-freedom. The system initially used 1500 control cycles per 
second and fixed point arithmetic programmed in assembly language. A sensitivity to input 
noise was discovered but these problems were overcome and development of a six degree-of- 
freedom digital controller began. 
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In 1984 a six deg rec-of- freedom digital controller was completed (Ref. 15). The system 
continued to use the digital phase-advance controller, with minor changes from the 1981 
algorithm in the first and third blocks. These changes were: 


f 1 "I 

_ t+ a s J; 


A T^ = f k _y k (P rcvious| y ; A T" = e k~ y k-i^ 


Ay 


where Ay = yj. - y k _, 


then y k - ( A ^ T ) < k + ( A + T ) ^k-i 


These equations can also be reduced to difference equations as: 


V' k = a 2 v Vi + a 3 a i £ k + a 4 a i f k-i 
V k = a 2^k-l"*" a 3 a l f, k'*' a/,aif, k-l 


(4.2a) 

(4.2b) 


Here aj= and a 2 = ^ A ,p which differ from equations 4.1a and 4.1b for the earlier systems. 

This form was believed to give superior performance for long sampling intervals (T»A). As 
extra control tasks placed increased demands on the control system, increased processing 
capability was necessary. This was provided by replacing the PDP-11/34 with a PDP-11/84. 
The extra control tasks included position sensor processing and output demand distribution 
related to high angle of attack operation (Ref. 1C). The control algorithm is in floating point 
assembly language and originally operated at 400 control cycles (all six degrees of freedom) per 
second. The controller now operates at 256 control cycles per second. 
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The error integrator used in the system is located in the forward path. The integrator in 


digital form is: 



The error integrator drives the steady-state error to zero. 


(4.3) 


4.1.4 NASA Langley, United States: 

In 1984 the NASA Langley Research Center 13-inch MSBS was converted to digital controls 
(Ref. 17). The controller closely followed the Southampton system, using a PDP-11/23 
computer to control five degrees-of-freedom (no roll control). With the same control loop 
configuration, the algorithm was modified slightly from the Southampton version to save time 
in execution (eliminated one floating point multiplication): 


f_> A _ “ TfHAS) y 


[ 1+ ii A 5]-* V' 


where y fc = + y k J and V' k - (1+4^) y k - ^ y k _, 


This allows the entire dual phase-advance transfer function to l>c rearranged as: 


c -* 


Tfl _ - e _ A 

|a 2 J [t(i+ a 


y -♦ Q+ nA S~] V' 



[l + nA S] -♦ V 
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This can be expressed as three equations applied in series as: 


e = a 2 c 

(4.4a) 

V' k = a 3 a,? k +a,* 4 C k . I +a I V , k . 1 

(4.4b) 

V k =a 1 a 3 V , k +a 1 a 4 V , k 1 +a,V k j 

(4.4c) 


where, a x = a 2 = a 3 = and a 4 = ^ 

The NASA controller uses floating point assembly language and a controller operating at 256 
cycles per second. 

4.1.5 NAL, Japan: 

The newest MSBS was commissioner! at the National Aerospace Laboratory, Tokyo, Japan 
in 1987, with digital controls used from the outset. Few details of the controller are available. 
However, the system appears to use some form of digital approximation of a classical P I D 
algorithm carried out on a microcomputer. 

r = _Kp + K d s+ 

Only three degrees-of-frecdom were controlled initially, but the system is designed and is being 
developed for full control of at least 5 and possibly 6 degrees of freedom (Ref. 18). 

A summary of digital controllers for MSBS wind tunnels is shown in table 2. 
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Table 2. Digital controllers For MSBS wind tunnels. 


Organization 

Date 

Degrees of 
Freedom 

Controller 

Type 

Oxford University 

1971 

3 

phase-advance 

MIT 

1976 

i 

phase-advance 

University of Southampton 

1981/84 

2/6 

phase-advance 

NASA Langley 

1984 

5 

phase-advance 

NAL 

1987/89 

3/5 or 6 

proportion- in tegral- derivative 


4.2 Other Digitally Controlled Magnetic Suspension Systems 

The first magnetic suspension system was originally developed for use as friction-free bearing 
for ultracentrifuge studies. Magnetic suspension systems are now being developed for 
transportation, magnetic bearings, and similar uses. It is worthwhile to review briefly the 
development of digital controllers for these uses since many of the problems and potential 
advantages are similar to those related to the use of digital cont rol for MSBS used with wind 
tunnels. 

4.2.1 Loughborough, England: 

A single dcgrec-of-freedom demonstration system has been developed at Loughborough 
University, England (Ref. 19). The digital controller algorithm approximates the output and 
input as quadratic curves^shown in figure 11. The controller is located in the forward path of 
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the circuit. The coefficients of these quadratic curves can then be solved if three points along 
the curve are known; the present, and two previous values. The quadratic curves are: 

input: c(t)= p + qt + rt 2 ( 

V(t)=P + Qt + /ft 2 


output: 




The coefficients of these quadratic curves can t>o solved in terms of the three points along 


the curves. The coefficients are: 


P = f k 


P = 


V k+1 (c 2 - 3C + 2 ) - 2 V k (C 2 - 20 + Vj^ <) 


q = 


(3v «k.,+‘k-,> P _ v k.,( 3 ■ 2 <) - " v k(l - 0 + v k-(» - 2 <> 
2T 


2T 


(4.6) 


(v *k-i+‘k.,> R = < V k + l - 2 V k + V k-,> 


2T 


2T 


To obtain the required control, the output + 1 of the controller is shifted forward an 
incremental time, This shift forward in time is called strike time and is designed to overcome 
calculation and system time lags. The controller provides a control command for a point in the 
near future. With these quadratic equations and an appropriate time shift forward for the 
output, the algorithm can represent several different controllers. The algorithm has the form: 


V k + 1= a ° f k + ^k-l + a -2 f k-2 + b ° V k + b -* V k-l 


(4.7) 


In equation 4,7 the coefficients are based on the quadratic curve fit coefficients obtained 
from equation 4.6 and the type of constants desired in the controller. If a dual phase-advance 
has the form j- = ( S ) ’ ^ en ^ ie coc ^ c ‘ ont,s c( l ua ^ on 4.7 arc: 


a o = (l+ 
a -i = (-2C 

a - 2 = ( 5 + 


K . 3nA , 
2 + T + 


c 2 

T 


2nAC 
+ T 


4nA ^2 4nA£ 2n 2 A 2 
-f" - k - T * T 2 


nA , C! , 2nAC , n?A_ 2 

T + 2 + T + x 2 
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4.2.2 Mitsui Engineering and Shipbuilding, Japan: 

For many years, researchers in Japan have studied the use of magnetic suspension for high- 
speed trains. The first known use of digital control techniques was with a magnetically 
suspended linear guide developed in 1984 by Nippon Telegraph and Telephone, Kanagawa ? as a 
technology demonstration (Ref. 20). Although the rate signals were derived from analog 
differentiators, the remainder of the control loop, including calculation of a coupling matrix, 
were carried out digitally. It appears that the digital hardware was custom built. 

A single degree-of-Creedom magnetic bearing has boon digitally controlled using a 
microcomputer by Mitsui Engineering and Shipbuilding, Okayama. Three approaches to the 
synthesis of the digital controller were tested (Ref. 21). 

The first approach used a digital simulation of a P I 1) controller. The digital controller 
uses the present and two previous position errors to determine the output command to the 
system. The rate prediction comes from a quadratic fit to the position error data. The values 
used for the proportional and derivative of the position error are at time 1.5T. 

T = K K p + K d 

Using a quadratic fit to the position error data as in the Loughborough system: 
f(t) = p+ ?t+rt 2 
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the position error and its derivative are calculated at time 1.5T. 


c(1.5T) = p + q (1.5T) + r (1.5T) 2 i(1.5T) = q + 2 r (1.5T) 2 


Also assuming the integral term is the sum of the |>osition error data over time, the P I D 
controller is: 



A second controller uses a P D controller using the same system as in equation 4.8 and 


letting the integral gain, Kj n = 0. 



The third method is a z-transformation of the P 1 D analog controller; where the P I D is 
represented as: 


V 

c 


= K 


Kp T Kj T 



_ K p S + K , S 2 + K 


in 


S 


Using a four-point central difference approximation for the derivatives of t, then: (Ref. 22) 


c k-2~ 8 f k-i+ 8 f k + r f k + 2 
12 T 


■■ - f k-2+ l6f k-.- 30f k+ 16f k+r c k+2 

f 12 T 2 


The first derivative of V is approximated by backward-difference where: 


V = 


V, -V 


k k-i 


31 



The P I D controller has the general form: 


V _ a 2 2 2 -f a t 2 -f a 0 + a.^' 1 + a_ 2 * 


T “ 


1 - * 


where: 


K P K d 

2_ * IT ‘ ITT 


_ 8 K p 16 K , 
ll ~ 12 + 12 T 


30 K , 

a ° “ T K in ‘ 12 T 


8 K p , 16 K d _ K p K d 

a -i w + ITT - 2 ~ 12 * ITT 


This is presented as a difference equation: 


V k = a 2 f k+2 + ai f k+i + a ° f k + a -‘ f k-i + a - 2 f k-2 + V k-i 


(4.9) 


The values of and fj {+2 are calculated by using the quadratic approximation of the position 
error as used in equation 4.5. 


f k + 2= 6 f k - 8 f k-. + 3 f k-2 


f k+ 


T 3 f k ' 3 f k-i + e k-2 


When these values are substituted back into equation 4.9, the output is expressed in terms of 


the inputs, e. , and 2 - 


(6a 2 +3a 1 +a 0 ) f^d- (-8a 2 -3a)d-a.i) j+ (3a 2 +a,+a_ 2 ) fj.. 2 + j 


(4.10) 



The next simplification is to let the initial value of the cont roller output, V 0 — 0, then for 


k=l, 2,3,4,. .. 

Vj= (6a2+3a x Ta 0 ) Cj-f (-8a 2 -3aj Ta^) c 0 + (3a 2 +a,+a_ 2 ) t„ { 

V 2 — (6a 2 +3a 1 +a 0 ) c 2 "h ( - 8a 2 -3aiTa_t) (3a 2 -|-a|-|-a_2) ^o“h ^1 

V 3 = (6a 2 +3a 1 -fa 0 ) e 3 + (-8a 2 -3a 1 -ha_ 1 ) c 2 + (3a 2 +a,+a. 2 ) * 1 + V 2 
V 4 = (6a 2 +3a 1 -fa 0 ) (-8a 2 -3a 1 -|-a_i) c 3 + (3a 2 Ta 1 +a_ 2 ) V 3 


The equations above can then be rewritten as: 


Vi = (5a 2 + 2a x - a. r a_ 2 ) “ (3a 2 + a j+a_ 2 ) c i .+ ( a 2 + a t+ a o+ a -i+ a - 2 ) £ c j 

K j=o 


(4.11) 


This is a P I D digital controller using only two position-error data points. 


4.2.3 Oak Ridge Gaseous Diffusion Plant, United States: 

The Oak Ridge gaseous diffusion plant in cooperation with the University of Virginia 
developed a digital magnetic bearing system (Ref. 23). The fundamental approach of the 
controller is to generate an estimate of the derivative of the suspended object by real-time curve 
fitting of the position data. This single degree-of-freedom cont roller uses a polynomial least 
squares fit with exponential weighting to estimate the derivative in a P I D controller. The idea 
of using exponential weighting is that the data furthest back in time from the present should 
have the least effect on the output. The form of the P I D control algorithm is: 



K p + K d 
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This controller assumes that the input to the controller is a polynomial curve of order n. 


c(t) — aQ -{- a t t *4* &2t 4* agt + - - . + a^t 


(4.12) 


The coefficients of this polynomial arc found by a least squares ‘best fit’ with a weighting factor. 
These coefficients of the polynomial are changed by incremental amounts, 6aj. The incremental 
changes Sa^ are determined by the order of the polynomial and the value of the weighting 
function, and remain constant. The calculated values for 6aj are: 


for (n=l) 6a 0 =l - W 2 6a,=(l - YV) 2 

for (n=2) tfao=l - W 3 «aj=|(l - W) 2 (J + W) 

^a 2 =I(l - W) 2 

for (n=3) <5a 0 =l - W 4 £a,=|(l - W) 2 (ll + 14VV + 11W 2 ) 

«a 2 =(l- W) 3 (l + W) «a 3 =i(l-W) 4 

These incremental changes can be calculated for an n*'* 1 order system. 

The algorithm is used in two forms. One calculates a present time output and the other 
calculates a predicted time output. 


For the present time output the controller algorithm is: 

a) Decide the polynomial order, n and the weighting factor, W 

b) Calculate the incremental changes in the polynomial coefficients, 6a^ 

c) Calculate present position error, = r - 

d) Calculate the change in the errors from predicted, Cj.- a 0 





c) Sum up position error for integral term, 

f) Apply incremental change to coefficients, aj^aj+^ajAc where i=0 to n 

g) Calculate present control output, 

V = K [Kp. k + K d », + 2 K dj , 2 + K in E<] 

h) Return to c) 

The predicted time output calculates the output from the controller at one time step 
forward using the present coefficients and shifting them forward. This is accomplished by 
substituting t=t+T into equation 4.12, where T is one time unit. Then: 

f'(t) = a 0 + a x (t+T) + a 2 (t+T) 2 + a 3 (t+T) 3 + . . . + a n (t+T) n 

Collecting the coefficients the predicted polynomial is: 

e ; (t) — a Q-f a jt -f a 2 t T • * • T a ^t 

where: a f 0 = a o + a i + a 2 4- a 3 + a 4 4- . . ■ 

a^ — a t 4- 2a 2 + 3a 3 + 4a 4 4* . . . 

a^ 2 = a 2 3a 3 -4 6a 4 -4 * * » 

a^ 3 — a 3 4" 4a 4 

a ^n = a n 

The predicted time output algorithm is: 

a) Decide the polynomial order, n and the weighting factor, W 

b) Calculate the incremental changes in the polynomial coefficients, 6aj 


(4.13) 


(4.14) 
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c) Calculate present position error, = r - 

d) Calculate the change in the errors form predicted, Ae= c^- a ; 0 

e) Sum up position error for integral term, 

f) Apply incremental change to coefficients, aj=a|+^ajAe where i=0 to n 

g) Shift the polynomial coefficients forward one time set. 

h) Calculate the predicted control output, 

v = K f K p+ a' 0 + K p a 0 + a', + 2 K , a' 2 + K () a, + 2 Kj a 2 + K^l (4.15) 

i) return to c) 

These two controllers allow the operator to select any order polynomial and any weighting 
factor for the control algorithm. The controller is located in the forward path. This system 
provides the most involved controller of all those discussed, 

4.2.4 UVa Electrical Engineering, United States: 

Magnetic suspensions are used at two locations at the University of Virginia and digital 
controllers are being developed for use with these systems. 

One group is the Electrical Engineering Department which is studying the use of magnetic 
bearings for a rotating shaft (Ref. 24), The system digitally controls the magnetic bearing 
through a microcomputer using assembly language. The magnetic bearing system uses a digital 
P D controller located in the feedback path where S = ^ The general form of a P D 
controller is: 

T = [ K p + V] 
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When the ^-transform of the derivative is substituted into the general form^the equation is: 


v k= K P 'k + t(v 'k-,) < 4 - 1C > 

This digital controller in equation 4.16 is very simple and provides adequate control of the shaft 
with the magnetic bearings. The simplicity of this program allows very high computational 
speeds. 

4.2.5 UVa Nuclear Engineering and Engineering Physics, United States: 

The other group at UVa is in the Department of Nuc lear Engineering and Engineering 
Physics. Magnetic suspension is used in this department for experimental studies of gravitation 
and general relativity (Ref. 25). The controller is a digital P I I) located in the forward path 
and analog filters. The digital P I D is of the form: 

V = K Kpf + K d e + K jn y f dt 

where the derivative is calculated using the first two terms of a Taylor series, where: 



The integral term is derived by using the Trapezoidal Rule, where: 
k 

] ( dt - ?( f k+ c k-i) 
k-i 
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This integral value is then summed up over the entire time. 


Ef - E f + J( f k+ c k-i) 


The proportional, derivative, and integral 


terms are then added to obtain the controller output. 


V = K 



+ 



e k-i 


+ 



(4.17) 


This digital controller is simple and provides adequate control of a suspended sphere. 


A summary of digital controllers for magnetic suspension systems 


is shown in table 3. 


Table 3. Digita 

controllers for magnetic suspension systems. 

Organization 

Date 

Controller Type 

Method 

Loughborough University 

1986 

phase-advance 

quadratic fit 

Mitsui Engineering and 
Shipbuilding 

1984 

P I D and P D 

quadratic fit 
difference equation 

Oak Ridge Gaseous Diffusion 
Plant 

1986 

p i n 

exponential weighting 
with polynomial fit 

University of Virginia, 

Dept, of Electrical Engineering 

1987 

P D 

^-transformation 

University of Virginia 
Dept, of Nuclear Engineering 
and Engineering Physics 

1989 

p 1 1) 

difference equation 
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5. DIGITAL SIMULATION 


5.1 Derivation of Equations for Simulation 

The MSBS plant described by equation 2.10 is for a single degree-of-freedom MSBS. This 
plant can be discretized by several different methods. One method is the Tustin’s 
transformation which is only an approximation of a conversion between the 5-domain and the z~ 
domain (Ref. 26). For the Tustin’s transformation: 


_ 2 (* 

- T (* + iy 


where T is the sampling time. 


As given in equation 2.10b the MSBS plant is: 


AX = 


^yav + JV0+ fc 5 ) f 

m L m I A II / 

c 3 . c 2 m.CN , C (cR . Kx. VA 
6 + * U +m ; + (m L 1,1 L m J 


R K x 
L in 


For zero input force disturbance, (f=0), the Tustin’s transformation of this equation is: 


AX _ a 0 [l+ z~' + z 2 + *~ 3 ] 

AV b 0 + bj z 1 + b 2 2 2 + b 3 z 3 


(5.1) 


-R. 

where a n = — f 
u m L 


b - a N 3 . (1 NVR . CN . 2 /C_R . Kx _ VA . Mx 
b ° “ ItJ + VtJ VL +rT V + Tl m T, rn TlmT I L m 
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For the case when f^O the Tustin’s transformed MSBS plant is: 


— a o[^d~ £ ] + z 2 -f z 3 ]AV + [ai+ a 2 z jh a 3 z 2 + a 4 z 3 ] f 
b 0 l*i z * + b 2 z 2 “h ^3 * 3 


(5.2) 


where the additional coefficients are: 


* l_ (mT, + n?r) 

-(H-A) 


( 3 11 

\m L 


+ 


2 

mT 


) 



mT/ 


Using equation 5/2 as the discretized MSBS plant, a simulation can be designed for use on a 
microcomputer. This simulation will allow design work for development and comparison of 
control algorithms. 


5.2 Simulation Program 

The simulation program is written in the BASIC language (Ref. 27). As with most 
microcomputer languages, BASIC allows great flexibility in the type of control algorithms that 
can be implemented on microcomputers for use as MSBS controllers. Because most MSBS 
systems use microcomputers to control the system, the BASIC language program can be used on 
an MSBS system or transformed to another computer language for use. 
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The simulation program includes the digital controllers discussed in chapter 4. The 
simulation allows the parameters of the MSBS plant to be changed easily and to observe the 
effects these variations have on the system performance. Also, the simulation program has the 
ability to vary the parameters in the controller and the type or method of control used for the 
MSBS. The program allows two types of step inputs to the system, a position input and a force 
input. 

A standard simulation run starts with a unit step position input at simulation time t=0. At 
simulation time t=5 seconds a 10-unit step force input is commanded. The simulation run then 
stops at t=20 seconds. The program has a graphical display of the suspended body position 
trajectory. This graphical display can be scaled to provide a detailed view of the trajectory. 

The program also calculates and displays, above the trajectories, certain design parameters that 
can be used to compare systems’ performances. These design parameters are gains, rise time, 
peak times, settling times, overshoots, time, and position. The complete listing of the program 
is given in Appendix A. 

5.3 Representative Magnetic Suspension and Balance System 

The choice of a representative MSBS plant for use in the simulation program is critical in 
order to determine how different controllers perform. This representative MSBS should exhibit 
the same dynamic characteristics as a real plant. These characteristics are determined by the 
location of the poles. As shown in equation 2.10, the pole locations are influenced by many 
parameters of the system. Many technical papers have addressed the problem of plant model 
verification with experimental results. The model described in equation 2.10 is more complex 
than most linearized models. Comparisons between experimental results and linearized models 
show that the dynamics of a magnetic suspension system are accurately described. 

To obtain the desired MSBS plant dynamics, three poles are needed with locations similar to 
those shown in figure 6. Based on the relative location of the poles for a real MSBS plant, a 
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suitable choice for the poles of our representative MSBS arc: P t c^ -10, P 2 — -1, and P 3 ~ 1. To 
realize these poles, the parameters of the system, as descril>ed in equation 2. 10, are: 


R = 1 

L = 0.1 

II 

X 

S4 

*— « 
O 

II 

uT 

o 

1 

II 

w 

ii 

E 

C = 0 




With these parameters the actual pole locations of the representative MSBS are Pj= -9.9899, 
P 2 = -1.0056, and P 3 = 0.9955. The plant transfer function therefore is: 


-AV + 10 fl + 0.1 S) f 

_ \ / 

S 3 + 10 S 2 - 0.9 S- 10 


(5.3) 


The actual choice of pole locations for this representative plant are not completely random. 
Recent work at the NASA Langley 13-inch MSBS has been toward developing a mathematical 
model of the system. The early results show the actual system has pole locations similar to 
those chosen for the representative plant. Also several reports have shown that the linear 
approximations give a good representation of the MSBS dynamics (Ref. 4). 
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6. COMPARISON OF CONTROLLERS BY SIMULATION PROGRAM 


The simulation program can be easily used to compare the responses of different control 
algorithms on a representative MSBS plant. As shown in chapter 4, there are several 
philosophies of how a digital controller is derived and the type of controller. Each method has 
advantages and disadvantages with the final decision based on the desired system performance. 

Studying digital controller algorithms is best carried out with a computer simulation 
program. Standard control systems analysis will not show the difference caused when deriving a 
digital controller. These differences are brought about because of approximations made when 
converting an analog controller to a digital controller. With the simulation program, the exact 
method of how the controller is executed can be programmed. The simulation allows the 
method of control to be changed or modified for comparison and development. The main 
purpose of the simulation is to study the different controllers to determine the advantages and 
disadvantages of a particular control system and compare several of their performance 
characteristics. 

6.1 Location of Controller 

The two primary uses for magnetic suspension systems are for large gap suspension and 
small gap suspension. The difference between large or small gap is based on the distance 
between the electromagnets and the suspended body. 

Large gap systems include those associated with wind tunnels. These systems require 
position input commands to change model position and orientation during wind tunnel tests. 

The wind tunnel system must also maintain position and orientation when loads are being 
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applied to the model. On the other hand, magnetic bearing suspension systems are small gap 
systems. Due to the small gaps, they seldom require a position input command, being mainly 
required to maintain a fixed position under applied loads. 

The large and small gap systems also have a difference in the power requirements. The 
current used in maintaining the suspension is several times greater in a large gap system than a 
small gap system. For example the NASA 13 inch MSBS requires approximately 20 amps in 
each coil to suspend a model; whereas, the Loughborough MSBS and other magnetic bearing 
systems use less than 1 amp in the coil. 

The different requirements and power levels for these two systems has produced two classes 
of controllers. Most wind tunnel suspension systems have used the phase advance controller 
located in the feedback path of the control circuit and an integral term located in the forward 
path. Typically these controllers have performed well to the position inputs and force inputs. 

Most magnetic bearing suspension systems use a P I D controller located in the forward 
path of the control circuit. This forward path P I D controller responds well to force inputs 
but poorly to step position inputs. 

The poor performance of the forward path controller to a step position input is caused by 
the lead compensation located in the forward path. Given a step position input, the initial 
derivative term of the controller is very large which causes a large first overshoot. The large 
overshoot is not a problem for a bearing system because position inputs are not expected. The 
bearing shaft would only momentarily touch the wall of the bearing and would quickly recover 
and continue to function properly. This large overshoot can be avoided by not allowing step 
inputs to the controller but rather limit the commands to ramp inputs. An advantage of having 
the controller in the forward path is to provide a quick response when compared to controllers 
located in the feedback path. 

An example of this large overshoot is shown in figure 12. The top graph is the position 
trajectory of a P I D controller with the P D part of the controller being located in the feedback 
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path. The lower graph is of the same P I D controller with the entire controller located in the 
forward path. Each of these controllers have identical gains and are subject to the standard 
simulation run. The only difference is the locat ion of the P and D parts of the controller. These 
gains are based on a 5% overshoot performance for a position input of a P I D controller with 
the P I) locat 'd in the feedback path. 



Time 

Figure 12. Position trajectories of P I D, (controller location). 
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Figure 12 shows the large overshoot produced by a step position input to the forward path 
controller as compared with the feedback controller. The response to the force inputs are almost 
identical with only minor differences caused by the integral term. The results of these responses 
are shown in table 4. The extremely large overshoot of 48% for the forward path controller is 
unacceptable. It is possible to design the forward path P I D controller having a 5% overshoot 
to a step position input but the controller then does not have sufficient stiffness to withstand 
large force inputs. By adjusting the overall gain pf the forward path P I D controller, a 
minimum First overshoot can be found. In table 4 the P I D “best is the best response to a unit 
step position input for the forward path P l D controller. 


Table 4. PIP controller location. (Position input) 


Controller 

Overall 

Gain 

Rise 
Time, s 

Peak 
Time, s 

Settling 
Time, s 

First 

Overshoot 

P I D, feedback 

362 

0.41 

0.84 

1.29 

1.050 

P I D, forward 

362 

0.17 

0.46 

1.15 

1.481 

P I D, “best” 

522 

0.13 

0.35 

1.64 

1.470 


This large overshoot is also present when using a dual phase-advance controller in the 
forward path. Figure 13 shows two position trajectories for dual phase-advance controllers. The 
top trajectory is for the controller located in the feedback path and in the bottom trajectory the 
controller is located in the forward path. Each controller is subject to a standard simulation 

run. 

Again the forward path has an unacceptable First overshoot. f l hese dual phase advance 
controllers are identical except for the location (>r the controller. The gains are based on a 5% 
overshoot performance for the feedback dual phase~ad vancc controller. The results of these 
responses for a stop position input are shown in table 5. The responses to the force input are 
nearly identical. 
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Table 5. D P A controller location. (Position input) 


Controller 

Overall 

Gain 



Settling 
Time, s 

First 

Overshoot 

D P A, feedback 

2080 

0.21 

0.-17 

0.93 

1.050 

D P A, forward 

2080 

0.07 

0.18 

0.63 

1.528 

D P A, “best” 

1602 

0.08 

0.22 

1.05 

1.487 


By adjusting the overall gain of the forward path dual phasc-advance controller, a minimum 
first overshoot can be found. In table 5, the dual phase advance “Iwst 1 ’ has the best response to 
a unit step input for the forward path controller. 

The rise, peak, and settling times are greatly improved by having the controller located in 
the forward path. However, these advantages are overshadowed by the unacceptable first 
overshoot. 

Using a controller in the forward path of a wind tunnel system could be dangerous. During 
a large overshoot the model could be lost from the view of the position sensors causing loss of 
model control. This is not to say that forward path controllers should never be used. However, 
care should be taken in the type of position inputs given to the controller. 

6.2 Comparison of Dual Phase Advance Controllers 

To compare the different algorithms of the digital phase-advance controllers, the constants 

within the controllers must be the same. Bach dual phase-advance controller is located in the 

feedback path and has an integrator added to the forward path to help improve performance by 

driving the steady-state error to zero. This integrator is based on equation 4.3. The integral 

gain is set at Kj =0.5 in all the algorithms. The controller time constant is also fixed at 

A=0.O1. The higli/low frequency gain is set to n=10. The only adjustable constant in each 

controller is the overall gain, K. With the constants being the same in each controller, the 

differences in performance of the digital dual phase-advance controllers can be compared. 
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There arc four dual phase-advance controllers that are compared using the simulation 
program. These arc the Tuslin’s D P A, Southampton, NASA, and Loughborough controllers. 
For the Tustin’s D P A algorithm, the Tuslin’s method is user! to discretize a dual phase- 
advance controller. The derivation of the Tustin’s D P A controller is shown in Appendix B as 
equation B-3. The Southampton controller is obtained from equations 4.2a and 4.2b, the NASA 
controller is from equations 4.4a, 4.4b, and 4.4c, and the Loughborough controller is from 
equation 4.7. 

Because the wind tunnel type controllers are concerned with position and force inputs, the 
performance analysis must include these inputs. A standard comparison run starts at the 
simulation time t=0 with a unit step position input. At simulation time t=5 seconds a 10 unit 
step force input is commanded. The computer program stops after 20 seconds of simulation 
time. The performance of the controller can be determined from these two input commands. 

6.2.1 5% Overshoot Per fo rmance: 

For dual phase-ad vance controllers, one design criterion for comparing the controllers is to 
adjust the overall gain for a first overshoot of 5% for a unit step position input. Figure 14 
shows the position trajectories for this 5% position input overshoot of each controller. 

The results of these trajectories arc shown in tables 6(a) and 6(b). Table 6(a) shows several 
performance parameters obtained from a position input. Table 6(b) shows the performance 
parameters obtained from a force input. 


49 



1.5 


C 





0 5 



Figure 14. Position trajectories of 0 P A contr 
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Tustin's D P A 



oilers, 5 % overshoot. 







Tabic 6(a). D P A, 5% overshoot. (Position input) 


1 rtuit; v 

Method of 
discretization 

tyaj, j-v » -t* 

Overall 

Gain 

^ ~ 

Rise 
Time, s 

Peak 
Time, s 

Settling 
Time, s 

First 

Overshoot 

Tustin’s D P A 

2784 

0.18 

0.36 

0.83 

1.050 

Southampton 

2665 

0.17 

0.34 

0.96 

1.050 

NASA 

2720 

0.17 

0.33 

0.95 

1.050 

Loughborough 

2080 

0.21 

0.47 

0.93 

1.050 


Table 6(b). D P A, 5% overshoot. (Force input) 


J. dUlt V 

Method of 
discretization 

Overall 

Gain 

, vrv 

Peak 
Time, s 

Settling 
Time, s 

Overshoot 

Final 

Position 

Tustin’s D P A 

2784 

0.24 

15.00 

1.083 

1.00 

Southampton 

2665 

0.23 1 

15.29 

1.090 

1.00 

NASA 

2720 

0.22 

11.34 

1.086 

1.00 

Loughborough 

2080 

0.27 

15.77 

1.102 

t.oo 


These results show that all the controllers perform well in controlling the system with li< in- 
difference in their performances. However the results for a position input show that the Tustin s 
D P A controller performs “best” because of its low settling time. The rise and peak times of 
the Tustin’s D P A controller are similar to those of the NASA and Southampton controllers. 

Table 6(b) shows the results from a force input. This is important because it shows the 
spring-like stiffness of the system which is caused by the controller. This stiffness is related to 
the overshoot caused by a force input. The Tustin’s D I* A and NASA controllers have nearly 
equal stiffness. Table 6(b) shows that the NASA method provides the “best” settling time from 

a force input to the controller. 

The integral gain, K jn has a major influence in the response to a force input. A high 
integral gain improves the response to force inputs by reducing the settling time. This high gain 

also increases instability. 
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6.2.2 Minimum First Overshoot, Performance: 


With increased overall gain for the dual phase-ad vance controller, the system response 
represents an overdamped system. The system performance is improved if the overall gain is 
increased so the first overshoot is minimum for a position input. Any increase in gain causes the 
second overshoot to be larger than the first overshoot. Figure 15 shows the position trajectories 
of each system based on this minimum first overshoot gain value. 

Tables 7(a) and 7(b) show the different controllers’ performances based on the minimum 
first overshoot system performance. 


Ta ble 7(a). D P A, minimum first overshoot. (Position input) 


A •A”/* 

Method of 
discretization 

Overall 

Gain 

Rise 
Time, s 

Peak 
Time, s 

Settling 
Time, s 

* 

First 

Overshoot 

Tustin’s D P A 

3202 

0.18 

0.35 

0.76 

1.005 

Southampton 

2939 

0.17 

0.33 

0.88 

1.014 

NASA 

2995 

0.17 

0.32 

1.04 

1.014 

Loughborough 

3115 

0.26 

0.54 

0.59 

1.013 


Table 7(b). 

D P A, minimum first overshoot. (Force in 

put) 

Method of 
discretization 

Overall 

Gain 

Peak 
Time, s 

Settling 
Time, s 

Overshoot 

Final 

Position 

Tustin’s D P A 

3202 

0.22 

14.80 

1.071 

1.00 

Southampton 

2939 

0.21 

14.89 

1.081 

1.00 

NASA 

2995 

0.21 

11.21 

1.077 

1.00 

Loughborough 

3115 

0.19 

14.84 

1.062 

1.00 


Tables 7(a) and 7(b) show that operating the system with minimum first overshoot 
improves the performance when compared to the 5% overshoot system shown in tables 6(a) and 
6(b). The minimum first overshoot controllers have better rise, peak, and settling times plus an 


increase in the stiffness of the system. 









From table 7(a), using the settling time as the performance criteria, the Loughborough 
controller is the “best”, except for its rise time and peak times. If choosing a controller based on 
the minimum first overshoot system performance, the Loughborough controller is preferred 
because of its low settling time. 

In table 7(b) the response of a force input is given which shows that the Loughborough 
controller provides the “best” stiffness. The NASA controller provides the “best” settling time. 

6.2.3 Execution Times: 

One of the seldom mentioned design criteria for digital controllers for MSBS is execution 
time. Execution time is extremely important in providing a good controller. Execution time is 
based on the number of calculations each controller must make in order to operate. Execution 
time is highly dependent upon the hardware and software that the controller uses. The faster 
controllers having a short execution time arc preferred if they can provide adequate control of 
the system. Table 8 shows a representative time required to complete 25000 cycles of the 
controllers using the simulation program. Each of these dual phase-advance controllers have 
nearly equal execution times. The best controllers in terms of execution time are the Tustin’s 
and Loughborough controllers. 


Table 8. D P A relative execution times. 
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The difference in the “best” controllers for a 5% overshoot system and a minimum first 
overshoot system shows how important it is to choose a design criteria which best suits the 
desired system performance. Because there are many possible uses for the controllers, one should 
base the choice of controller on the expected use and desired performances of the system. 

6.3 Comparison of Proportional Integral Derivative Controller 

As shown earlier, controllers located in the forward path of an MSBS system have a large 
first overshoot to step- position input. Comparison of controllers in their forward path is useful 
even though they would perform better if located in the feedback path. Most of the P I D 
controllers used with magnetic bearing systems are located in the forward path. The comparison 
of the forward path controllers is useful because the response to force inputs and the stiffness of 
the system will be the same regardlcsss of the controllers’ location. These force responses and 
stiffness can then be compared with other controllers. 

To compare the P I I) controllers in the forward path, a set of design criteria must be 
established. To compare the P I D controllers, each controller must have the same gains within 
the controller and adjust only the overall gain of the controller. The constants chosen are 
Kp=l, 1^=0. 4, and Kj n =0.5. The value of the integral gain is the same as the value used in 
the phase advance controller comparison. The proportional and derivative gains were obtained 
from an analysis of the UVa controller described in equation 4.16. The proportional and 
derivative gains are based on the l>cst possible response of this P I D controller. As shown 
earlier, a P I D controller in the forward path does not have an acceptable first overshoot 
response to a unit step position input. 

There are six P I D controllers which are compared using the simulation program. These 
P I D controllers will be referred to as Tustin’s P I D, equation 4.8, equation 4.10, equation 
4.11, equation 4.16 and equation T17. The Tustin’s P I D controller is derived in Appendix B 
as equation B-5. The controller described by equations 4.8, 4.10, and 4.11 are all derivations 
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from Mitsui Engineering and Shipbuilding. The fifth controller, equation 4.16, is from the UVa 
Electrical Engineering Department and equation 4.17 from theTJVa Nuclear Engineering and 
Engineering Physics Department. 

One controller described in chapter 4 that is not used in this comparison is the Oak Ridge 
controller described by equations 4.13 and 4.15. The described controller algorithm will not 
control the simulated MSBS plant. The published documentation of this control algorithm is 
thought to contain an error and further clarification is being sought (Ref. 23). 

6.3.1 Minimum First Overshoot Performance: 

The design criteria for comparing the P I D controllers are to minimize the first overshoot 
for a unit step-position input and to have the highest possible overall gain. In each case there is 
a unique gain which provides a minimum first overshoot. 

The simulation run is the same as used with the dual phase advance controllers. At 
simulation time t=0 a unit stcj>-position input is commanded. Following this at time t=5 
seconds a step-force input of 10 units is given. The simulation stops after 20 seconds. 

The resulting position trajectories for each controller are shown in figure 16. The 
performance characteristics of these P I D controllers are presented in tables 9(a) and 9(b). 


Table 9(a). PIP, minimum first overshoot. (Position input) 


Method of 
discretization 

Overall 

Gain 

Rise 
Time, s 

Peak 
Time, s 

Settling 
Time, s 

First 

Overshoot 

Tustin’s P I D 

390 

0.16 

0.43 

1.38 

1.569 

Equation 4,8 

522 

0.13 

0.35 

1.64 

1.470 

Equation 4.10 

359 

0.17 

0.46 

1.45 

1.594 

Equation 4.11 

360 

0.17 

0.46 

1.45 

1.594 

Equation 4.16 

361 

0.17 

0.46 

2.02 

1.598 

Equation 4.17 

359 

0.17 

0.46 

1.47 

1.598 
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Table 9(b). PIP, minimum first overshoot, (Force input) 


Method of 
discretization 

Overall 

Gain 

Peak 
Time, s 

Settling 
Time, s 

Overshoot 

Final 

Position 

Tustin’s P I D 

396 


15.71 

1.506 

1.00 

Equation 4.8 

522 

0.55 

15.02 

1.346 

1.00 

Equation 4.10 

359 


15.57 

1.567 

1.00 

Equation 4.11 

360 

■ 

15.43 

1.565 

1.00 

Equation 4.16 

361 

0.62 

15.39 

1.568 

1.00 

Equation 4.17 

359 

0.63 

15.47 

1.573 

1.00 


From figure 16, and tables 0(a) and 9(b) the “best” FID cohtroller is described by 
equation 4.8. This controller has the fastest rise and peak times and the lowest overshoot for a 
position input. The Tustin’s P I D has the “best” settling time. For a force input, the equation 
4.8 controller has the largest stiffness as shown by the low overshoot value from a force input as 
shown in table 9(b). The equation 4.8 controller also has the “best” settling time from a force 
input. 

The controllers of equations 4.11) and 4.1 1 are nearly identical in response because equation 
4.11 is a derivation of equation 4. 10. Tin- method used in deriving the algorithm for equations 
4.16 and 4.17 is simple and provides a similar response to the complex algorithms of equations 
4.10 and 4.11. 
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6.3,2 Execution Times: 


The execution times for the P I D controllers were determined by the same method as used 
for the execution times of the dual phase advance controllers. Table 10 shows the execution 
time required to complete 25000 cycles of the controller. The table shows that each P I D has 
nearly identical execution times. 


Table 10. P I D relative execution times. 


Method of 
discretization 

Execution 
Time, s 

Tustin’s P I D 

35 

Equation 4.8 

35 

Equation 4.10 

34 

Equation 4.11 

34 

Equation 4.16 

35 

Equation 4.17 

35 
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7. CONCLUSIONS AND RECOMMENDATIONS 


7.1 Best Overall Controller 

The choice of the “best” overall controller is completely dependent on the desired system 
performance and intended use. The best location for the controller is the feedback path. The 
advantages of a forward path controller in reducing the rise, peak, settling, and execution time, 
do not overcome the inability to adequately control the system for a step-position input. The 
dual phase advance controllers provide superior performance in controlling the representative 
MSBS system when compared to the P I D controller. The choice of the feedback dual phase- 
advance controller as the “best” is based on the controller's suitability for a large gap MSBS 
system. The dual phase-advance controllers provide better stiffness than the P I D controllers. 

The “best” of the forward path P I D controllers is the Mitsui Engineering and 
Shipbuilding, equation 4.8. This controller is derived using a quadratic fit. This type of 
quadratic fit transformation also produced the “best” overall feedback dual phase-advance 
controller from Loughborough, equation 4.7. These quadratic fit controllers are simple to derive 
when compared to some of the other controllers. A feature of the quadratic fit is the selection of 
strike time, ^ which influences the response of the system. This strike time can be chosen to 
optimize a given system’s performance. 
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Nearly all the controllers have the same basic generic equation as shown below. 


V k = ^ C k + 2 + a l f k + l + a 0 f k + H -* <k-l + a -2 f k-2 + b -' V k-1 + b ‘2 V k-2 


The only difference is the method used to determine the coefficients and the coefficient values. 
The quadratic fit controllers provide both good control and a simple method of deriving the 
controller coefficients. The values of e k+2 and e k+J are controller input values which are future 
values that have not occurcd. These values are predictive by the quadratic fit controllers. 


7.2 Future Methods of Control 

With the development of modern control theories, the application of state-space type 
controls to an MSBS is likely to be an extension for future controllers. As shown earlier, the 
linearized mathematical model of an MSBS is both observable and controllable. This allows the 
selection of any desired system performance by the pole placement methods. These pole 
locations are only limited to the ability of the power supply. Another advantage of a state-space 
controller is in the simplicity of implementing the controller algorithm on computers. As with 
digital simulation of analog controllers, the possibilities of state-space controllers are also 
unlimited. 

One of the requirements for MSBS systems is the feedback signal to obtain stability. This 
feedback signal is usually body position, which is used to determine a velocity /derivative control 
signal. The idea of using acceleration feedback which can be integrated to obtain velocity and 
position is possible. The instrumentation to produce this feedback must be adaptable to strong 
magnetic fields. ONERA, in 19C8, suspended a model with a telemetry package that included 
four accelerometers (Ref. 28). The response times of the controllers will improve using 
acceleration feedback. Work is presently underway to study the use of acceleration feedback in 
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an MSBS system. The final goal would he the development of an internal rate gyro to obtain 
all the position information of the model. 

AH the controllers discussed in this report are linear controllers which do provide adequate 
control of an MSBS system. An improvement in system performance can be obtained by the 
development of nonlinear or adaptive controllers. These controllers will be more complex to 
develop and program. Presently, some digital controllers do have nonlinear controls which limit 
the output command to the power supply so as not to exceed its capabilities. The need for 
nonlinear controllers is evident in the wind tunnel because of the large changes in forces or body 
orientation during a run. With the present controllers, a standard wind tunnel run requires the 
operator to change the controller gains when the forces on the body change. These gain changes 
are referred to as gain scheduling and have been used at the University of Southampton in 
obtaining high angle-of-attack suspension. (Ref. 16) 

7.3 Effects that any Approximations m ay have on Results 

Several approximations are made in the derivation of the governing equation for an MSBS 
system. These approximations are considered reasonable simplifications to the nonlinear 
equations of a true MSBS. Several reports have shown that the linear approximation gives a 
good representation of the MSBS dynamics (Ref. 4). These approximations apply well to the 
magnetic bearing systems and to the wind tunnel systems while operating at their equilibrium 
points. The equations do not adequately represent the dynamics during large position changes 
away from the equilibrium point. In practice, the controllers which are designed using the linear 
MSBS plant also adequately control the system during large position changes from equilibrium. 

For any MSBS, the choice of a controller is extremely important because the controller will 
directly determine the performance of the system. However, the most important choice for any 
MSBS system is the available power supply. An ideal controller can have output commands 
that are beyond the capabilities of the power supply. It is possible to operate a system where 


62 


the power supply capabilities are low. Great care is required in the type of commands or loads 
applied to such a system. The limitations of the power supply arc not usually a problem with 
magnetic bearings because of the low currents used. For the large gap MSBS, as in a wind 
tunnel, the power supply limitations are a continuing concern. The limitations of the power 
supply used to provide the required currents to the suspension coils have not been covered 
extensively. This could allow a controller to be chosen as the “best” which requires more power 
than is available. A designer should be constantly aware in the choice of the best controller. 

7.4 Applications to Multi-Degree of Freedom System 

In a multi-degree-of-freedom system, several controllers must act together to maintain 
stability. For a multi-degree-of-frecdom system, the relation of the magnetic forces to body 
position are highly coupled and largely dependent on the arrangement of the coils. Decoupling 
of this relation into the required degrees-of-freedom is required for control. This decoupling is 
presently done for all MSBS associated with wind tunnels with good results in controlling a 
specific degree-of freedom. There is a slight coupling between some degrees-of-freedom; however, 
this quickly dies out in a few computational cycles. 
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APPENDIX A 


Program Listing 

The simulation program is written in Microsoft Quick Basic, Version 4.5. Below is a block 
diagram of the controller and MSBS plant used in this program. 




This program displayed and saved the position trajectories of a simulation run. Below is a 
printout of' the program’s display. 


KSBS Sinulation 

11:44:54 Tustin’s Method, Dual Phase Advance 

K=-27B4 


Tine = 16.086 
Rise Tine = 0.188 
Peak Tine 1 = 5.248 
Settling Tine 1 = 8.838 
Overshoot 1 = 1.883 


Position = 1.886 

Peak Tine 2 = 5.240 
Settling Tine 2 = 0.888 
Overshoot 2 = 1.083 



(Program Listing) 

CLS 

CLEAR 

’Saved as MSBSSIM.BAS 
’OPEN ”B:filename” FOR OUTPUT AS #1 
’Sampling Time 
T = .01 

’The MSBS plant variables 
Kx = -1 
Kc = -.1 
Ki = .1 
m = 1 
R = 1 
L = .1 
C = 0 

’The MSBS plant coefficients 
aO = -Ki / m / L 
al=R/m/L + 2/ m/T 
a2 = 3*R/m/L + 2/ m/ T 
a3 = 3*R/m/L-2/m/T 
a4 = R/ m/ L- 2/ m/ T 

bO = (2/T)“3+(2/T)‘2*(R/L+C/m)+2/T*(C*R/m/L+Kx/in-Ki*Kc/L/m)+R*Kx/L/m 
bl =-3*(2/T)‘3-(2/T)‘2*(R/L+C/m)+2/T*(C*R/m/L+Kx/m-Ki*Kc/L/m)+3*R*Kx/L/m 
b2 = 3*(2/T)'3-(2/T)"2*(R/L+C/m)-2/T*(C*R/m/L+ Kx/rn- Ki*Kc/L/m)+3*R*Kx/L/m 
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b3 = -(2/T)\3+(2/T)~2*(R/LHXym)-2/T^C*R/m/L+Kx/m-Ki*Kc/L/m)+R*Kx/L/m 
’Screen layout 
tmax = 20 
tmin = 0 


XMAX = 2 
xmin — 0 


SCREEN 9 
COLOR 14, 1 

VIEW (40, 125)-(620, 320), 9 

WINDOW (tmin-.01*tmax, xmin-.02*XMAX)-(tmax-f .OUtrnax, XMAX+.02+XMAX) 
’Borders 

LINE (tmin, xmin)-(tmin, XMAX), 14 
LINE (tmin, xmin)-(tmax, xmin), 14 
LINE (tmin, XMAX)-(tmax, XMAX), 14 
LINE (tmax, xmix)-(tmax, XMAX), 14 
’Horizontal lines 

LINE (tmin, .25 * XMAX)-(tmax, .25 * XMAX), 8, , &IIFF00 
LINE (tmin, .5 * XMAX)-(tmax, .5 * XMAX), 11, , &IIFF00 
LINE (tmin, .75 * XMAX)-(tmax, .75 * XMAX), 8, , &IIFF00 
’Vertical lines 

LINE (tmax * .1, xmin)-(tmax * .1, XMAX), 8, , &IIFF00 
LINE (tmax * .2, xmin)-(tmax * .2, XMAX), 8, , &JTFF00 
LINE (tmax * .3, xmin)-(tmax * .3, XMAX), 8, , &HFF00 
LINE (tmax * .4, xmin)-(tmax * .4, XMAX), 8, , &1IFF00 
LINE (tmax * .5, xmin)-(tmax * .5, XMAX), 8, , &IIFF00 
LINE (tmax * . 6 , xmin)-(tmax * . 6 , XMAX), 8, , fcHFFOO 
LINE (tmax * .7, xmin)-(tmax * .7, XMAX), 8, , &HFF00 
LINE (tmax * .8, xmin)-(tmax * .8, XMAX), 8, , &TIFF00 
LINE (tmax * .9, xmin)-(tmax * .9, XMAX), 8, , &IIFF00 
’Label 

LOCATE 12, 2: PRINT ”P” 

LOCATE 13, 2: PRINT ”o” 

LOCATE 14, 2: PRINT ”s” 

LOCATE 15, 2: PRINT ”i” 

LOCATE 16, 2: PRINT ”t” 

LOCATE 17, 2: PRINT ”i” 

LOCATE 18, 2: PRINT ” 0 " 

LOCATE 19, 2: PRINT ”n” 

LOCATE 2, 30: PRINT ”MSBS Simulation” 

’Input step of position 
ref = 1 


LOCATE 3, 2: PRINT TIMES 


Total = Total 
’GOSUB 100 
’GOSUB 250 
’GOSUB 300 
’GOSUB 450 
’GOSUB 475 
’GOSUB 525 
’GOSUB 550 
’GOSUB 600 


+ T 

’Tustin’s D P A, equation (B-3) feedback path 
’NASA D P A, equations (4.4a), (4.4b), and (4.4c) feedback path 
’Southampton D P A, equations (4.2c), and (4. 2d) feedback path 
’Loughborough D P A, equation (4.7) feedback path 
’Loughborough D P A, equation (4.7) Feedback path 
’UVa P D, equation (4.16) feedback path 
’UVa P D, equation (4.16) forward path 
’Japan P I D, equation (4.8) forward path 
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’GOSUB 650 ’Japan P D, and I, equation (4.8) feedback path 
’GOSUB 750 ’Japan P I D, equation (4. 11) forward path 

’GOSUB 800 ’Japan P I D, equation (4.10) forward path 

GOSUB 850 ’UVa P I D, equation (4.17) forward path 

’GOSUB 900 ’Oak Ridge P I D, equation (4.13) forward path, (Not Working) 

’GOSUB 1000 ’Oak Ridge P I D, equation (4.15) forward path, (Not Working) 

’GOSUB 1100 ’Tustin’s P I D, equation (B-5) forward path 
’Total Error Sum 

SUError = ABS(Xpl - X) / ref + SUError 
’Max. Overshoot and Peak Time for Position Input 
IF X > MAXX THEN MAXX = X 
IF X = MAXX THEN PTIME = Total 
’Max. Overshoot and Peak Time for Force Input 
IF Total > 5 AND X > MAX2 THEN MAX2 = X 
IF X = MAX2 THEN PT1ME2 = Total 
’Rise Time 

IF X <= (.1 * ref) THEN RT1 = Total 

IF jj = 1 THEN GOTO 98 

IF X >= (.9 * ref) THEN jj = 1 

IF X >= (.9 * ref) AND jj s 1 THEN RT2 = Total 

RISE = RT2 - RT1 

LOCATE 6,15: PRINT ”Risc Time =": LOCATE 6,27: PRINT USING ”##«###”; RISE 
98 LOCATE 9,13: PRINT ” Overshoot 1 =”: LOCATE 9,27: PRINT USING ”##.###”; MAXX 
LOCATE 7,13: PRINT "Peak Time 1 =”: LOCATE 7,27: PRINT USING "##.###”; PTIME 
LOCATE 7,43: PRINT "Peak Time 2 =”: LOCATE 7,57: PRINT USING ”##.##”; PTIME2 
LOCATE 5,20: PRINT "Time =": LOCATE 5,26: PRINT USING " ##.###”; Total 
LOCATE 5,46: PRINT "Position =”: LOCATE 5,57: PRINT USING ”##.###”; X 
LOCATE 9,43: PRINT "Overshoot 2 LOCATE 9,57: PRINT USING "##.###”; MAX2 
’Position Input Settling Time 

p = .001 

IF jjj = 1 THEN GOTO 59 

IF (ABS(Xpl-X)<p*X) AND (ABS(Xp2-X)<p*X) AND (ABS(Xp3-X)<p*X) AND (ABS(Xp4- 
X)<p*X) AND (ABS(Xp5-X)<p*X) AND (ABS(XP6-X)<p*X) AND (ABS(Xp7-X)<p*X) AND 
(ABS(Xp8-X)<p*X) THEN SETTIME = Total 
IF SETTIME = Total THEN jjj = 1 
LOCATE 8,9: PRINT"SettlingTime 1 = ": 

LOCATE 8,27: PRINT USING " ##.### ”;SETTIME 

59 ’Force Input Settling Time 
pp = .0005 

IF jjjj = 1 THEN GOTO 70 

IF fd > 1 THEN GOTO 60 ELSE GOTO 70 

60 IF ABS(X - ref) / ref < pp AND ABS(Xpl - ref) / ref < pp AND ABS(Xp2 - ref) / ref < pp 
AND ABS(Xp3 - ref) / ref < pp AND ABS(Xp4 - ref) / ref < pp AND ABS(Xp5 - ref) / ref < 
pp AND ABS(XP6 - ref) / ref < pp AND ABS(Xp7 - ref) / ref < pp AND ABS(Xp8 - ref) / ref 
< pp THEN SETTIME2 = Total 

IF SETTIME2 = Total AND Total > 6 THEN jjjj = 1 
LOCATE 8,39: PRINT "Settling Time 2 =”: 

LOCATE 8,57: PRINT USING "##.### SETTTME2 
’Shift the variables back in time 
70 fdp3 = fdp2 
fdp2 = fdpl 
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fdpl = fd 
Xp8 = Xp7 
Xp7 = XP6 
XP6 = Xp5 
Xp5 = Xp4 
Xp4 = Xp3 
Xp3 = Xp2 
Xp2 = Xpl 
Xpl = X 
Ep3 = Ep2 
Ep2 = Epl 
Epl = E 
Vp3 = Vp2 
Vp2 = Vpl 
Vpl = V 

PSET (Total, X), 15: 

’PRINT #1, USING ” ###.###”; Total; X 
’Input step of force 
IF Total > 5 THEN fd = 10 

IF Total > tmax AND Total < tmax + T THEN GOTO 80 ELSE GOTO 52 
80 LOCATE 24, 37: PRINT ’’Time” 

BEEP 

LOCATE 3, 2: PRINT TIMES 
CLOSE 
88 END 

’Subroutines 

100 ’Tustin’s D P A, equation (B-3), feedback path, plus error integrator 
IF first = 1 GOTO 110 
K = -2784 

LOCATE 3, 21: PRINT "Tustin’s D P A, equation (B-3)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kin = .5 

A = .01 

n = 10 

cO = (T*T + 4*n*A*T + 4*n*n*A*A)/(T*T + 4*A*T + 4*A*A) 
cl = (2*T*T-8*n*n*A*A)/(T*T + 4*A*T + 4*A*A) 
c2 = (T*T-4*n*A*T + 4*n*n*A*A)/(T*T + 4*A*T + 4*A*A) 
c3 = (2*T*T-8*A*A)/(T*T + 4*A*T + 4*A*A) 
c4 = (T*T-4*A*T + 4*A*A)/(T*T + 4*A*T + 4*A*A) 
dl = Kin * T 
110 E = ref * gain - G 
Etotal = E + Etotal 
Z = E + dl * Etotal 
V = K * Z 

X = (a0*(V-pVpl+Vp2+Vp3)+al*fd+a2*fdpl-t-H3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
G = c0 * X + cl * Xpl -I- c2 * Xp2 - c3 * Gpl - c4 * Gp2 
Gp2 = Gpl 
Gpl = G 
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first — 1 
RETURN 

250 ’NASA D P A, equations (4.1a), (4.4b), and (4.4c), feedback path, plus error integrator 
IF first - 1 GOTO 260 
K = -2777 

LOCATE 3, 26: PRINT "NASA D P A, equations (4.4a), (4.4b), and (4.4c)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain — 1 

Kin — .8 

A = .01 

n = 10 

cl = A / (T + A) 
c2 = (T * T) / (A * A) 
c3 = l + n*A/T 
c4 = -n * A / T 
dl = Kin * T 
d2 = cl * c3 * c2 
d3 = cl * c4 * c2 
d4 = cl * c3 
d5 = cl * c4 

260 E = ref * gain - G 
Etotal = E + Etotal 
Z = E + dl * Etotal 
V = K ♦ Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl-fa3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2- b3*X P 3)/bO 

UU = d2 * X + d3 * Xpl + cl * UUpl 

G = d4 * UU + d5 * UUpl + cl * Gpl 

UUpl = UU 

Gpl = G 

first — 1 

RETURN 

300 ’Southampton D P A, equations (4.2c) and (4. 2d), feedback path, plus error integrator 
IF first = 1 GOTO 310 
K = -2665 

LOCATE 3, 22: PRINT "Southampton 1) P A, equation (4.2c), and (4. 2d)” 

LOCATE 4, 35: PRINT ”K="; K 

BEEP 

gain = 1 

Kin = .5 

A = .01 

n = 10 

cl = T / (A + T) 
c2 = A / (A + T) 
c3 = (T + n * A) / T 
c4 = -n * A / T 
dl = Kin * T 
d2 = c3 * cl 
d3 = c4 * cl 
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310 E = ref * gain - G 
Etotal = E + Etotal 
Z = E + dl * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2-|-a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
UU = c2 * UUpl + d2 * X + d3 * Xpl 
G = c2 * Gpl + d2 * UU + d3 * UUpl 
UUpl = UU 
Gpl = G 

first = 1 
RETURN 

450 ’Loughborough D P A, equation (4.7), feedback path, plus error integrator 
IF first = 1 GOTO 460 
K = -2080 

LOCATE 3, 19: PRINT ” Loughborough D P A, equation (4.7)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kin = .5 

zeta =1.5 

A = .01 

n = K) 

aL0=(l+3*zeta/2+3*n*A/T-|-zeta‘2/2+2*n*A*zota/T-Hi*n*A*A/T/T)/(l-f3*A/T+A*A/T/T) 
aLl=(-2*zeta-4*n*A/T-zeta'2-4*n*A*zela/T-2*n*n*A*A/T/T)/(l-|-3*A/T+A*A/T/T) 
aL2=(zeta/2-fn*A/T-t-zeta'2/2+2*n*A*zcta/T+n*n*A*A/T/T)/(l+3*A/T-(-A*A/T/T) 
bL0 = (4*A/T + 2*A*A/T/T)/(l+3*A/T + A*A/T/T) 
bLl=-l*(A/T + A*A/T/T)/(l + 3*A/T + A*A/T/T) 
dl = Kin * T 
460 E = ref * gain - G 
Etotal = E + Etotal 
Z = E + dl * Etotal 

V = K + Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a4*rdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
G = aLO * X + aLl * Xpl + aL2 * Xp2 + bLO * Gp + 1»L1 * Gpl 
Gpl = Gp 
Gp = G 

first = 1 
RETURN 

475 ’Loughborough D P A, equation (4.7), forward path, plus error integrator 
IF first = 1 GOTO 460 
K = -3115 

LOCATE 3, 19: PRINT "Loughborough D P A, equation (4.7)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kin = .5 

zeta =1.5 

A = .01 

n = 10 
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aL0=(l+3*zeta/2+3*n*A/T+zcta"2/2+2*n*A*zcta/T+n*n*A*A/T/T)/(l+3*A/T+A*A/T/T) 
aLl=(-2*zeta-4*n*A/T-zeta‘2-4*n*A*zeta/T-2*n*n*A*A/T/T)/(l+3*A/T+A*A/T/T) 
aL2=(zeta/2+n*A/T+zeta'2/2+2*n*A*zeta/T+n*n*A*A/T/T)/(l+3*A/T+A*A/T/T) 
bLO = (4*A/T + 2*A*A/T/T)/(1+3*A/T + A*A/T/T) 
bLl=-l*(A/T + A*A/T/T)/(l+3*A/T + A*A/T/T) 
dl = Kin * T 
495 E = ref + gain - G 
Etotal = E + Etotal 

ZZ = aLO * E + aLl + aL2 * Ep2 + bLO * ZZp + bLI * ZZpl 
Z = ZZ + dl * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 

ZZpl = ZZp 

ZZp = ZZ 

G = X 

first = 1 

RETURN 

525 ’UVa P D, equation (4.16), feedback path, plus error integrator 
IF first = 1 GOTO 535 
K = -405 

LOCATE 4, 35: PRINT ”K=”; K 

LOCATE 3, 20: PRINT "UVa P D, equation (4.16)” 

BEEP 
gain = 1 
Kp = 1 
Kd = .4 
Kin = .5 
cl = Kd / T 
c2 = Kp + cl 
dl = Kin * T 
535 E = ref * gain - G 
Etotal = E + Etotal 
Z — E -f dl * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)yb0 
G = c2 * X - cl * Xpl 
first = 1 
RETURN 

550 ’UVa P D, equation (4.16), forward path, plus error integrator 
IF first = 1 GOTO 560 
K = -361 

LOCATE 3, 22: PRINT ”UVa P D, equation (4.16)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain — 1 

Kp = 1 

Kd = A 

Kin = .5 

cl = Kd / T 
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c2 = Kp + cl 
dl = Kin * T 
5C0 E — ref * gain - G 
Etotal = E + Etotal 
ZZ = dl * Etotal 
Z = c2 * E - cl * Epl + ZZ 

V = K * Z 

X=(aO*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a'l*fdp3-bUXpl-b2*Xp2-b3*Xp3)/bO 
G = X 

first = 1 
RETURN 

600 ’Japan P I D, equation (4. 8), forward path 
IF first = 1 GOTO 610 
K = -522 

LOCATE 3, 15: PRINT "Japan P I D, equation (4.8)” 

LOCATE 4, 35: PRINT "K=”; K 

BEEP 

gain = 1 

Kp = 1 

Kd = .4 

Kin = .5 

cl = 15 * Kp / 8 + 2 * Kd / T 

c2 = 42 * Kp / 8 + 5 * Kd / T 

c3 = 35 * Kp / 8 + 3 * Kd / T 

dl = Kin * T 

610 E = ref * gain - G 
Etotal = E + Etotal 

Z =s cl * Ep2 - c2 * Epl + c3 * E + dl * Etotal 

V = K * Z 

X=(aO*(V+Vpl+Vp2+Vp3)+al *fd+a2*fdpl+a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
G = X 
first = 1 
RETURN 

650 ’Japan P D and I, equation (4.8), feedback path, plus error integrator 
IF first = 1 GOTO 660 
K = -391 

LOCATE 3, 15: PRINT "Japan, P D, and I, equation (4.8)" 

LOCATE 4, 35: PRINT "K=”; K 

BEEP 

gain = 1 

Kp = 1 

Kd = .4 

Kin = .5 

cl = 15 * Kp / 8 + 2 * Kd / T 

c2 = 42 * Kp / 8 + 5 * Kd / T 

c3 = 35 * Kp / 8 + 3 * Kd / T 

dl = Kin * T 

660 E = ref * gain - G 
Etotal = E + Etotal 
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Z = E + dl * E total 

V = K * Z 

X=(aO*(V+Vpl+Vp2+Vp3)-fal*fd-fa2*fdpl -fa3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 

X total = X + Xtotal 

G = cl * Xp2 - c2 * Xpl + c3 * X 

first = 1 

RETURN 

750 ’Japan P I D, equation (4.11), forward path 
IF first = 1 GOTO 760 
K = -360 

LOCATE 3, 18: PRINT "Japan P I D, equation (4.11)" 

LOCATE 4, 35: PRINT ”K="; K 

BEEP 

gain = 1 

Kp = 1 

Kd = .4 

Kin = .5 

aJ2 = -Kp / 12 - Kd / (12 * T) 
aJl = 8 * Kp / 12 + 16 * Kd / (12 * T) 
aJO = Kin * T - 30 * Kd / (12 * T) 
aJpl = -8 * Kp / 12 + 16 * Kd / (12 * T) 
aJp2 = Kp / 12 - Kd / (12 * T) 
cl = 5 * aJ2 + 2 * aJl - aJpl - aJp2 
c2 = 3 * aJ2 + aJl + aJp2 
c3 — aj 2 -J- aj 1 -1- aJ 0 -1- aJ pi -|- aJ p2 
760 E = ref * gain - G 
Etotal = E + Etotal 
Z = cl * E - c2 * Epl + c3 * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
G = X 
first — 1 
RETURN 

800 ’Japan P I D, equation (4.10), forward path 
IF first = 1 GOTO 810 
K = -359 

LOCATE 3, 15: PRINT "Japan P I D, equation (4.10)” 

LOCATE 4, 35: PRINT "K="; K 

BEEP 

gain = 1 

Kp = 1 

Kd = .4 

Kin = .5 

aJ2 = -Kp / 12- Kd / (12 * T) 

aJl = 8 * Kp / 12 + 16 * Kd / (12 * T) 

aJO = Kin * T - 30 * Kd / (12 * T) 

aJpl = -8 * Kp / 12 + 16 * Kd / (12 * T) 

aJ P 2 = Kp / 12 - Kd / (12 * T) 

cl = 6 * aJ2 -(- 3 * aJl + aJO 
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c2 = -8 * aJ2 - 3 * aJl + aJpl 
c3 = 3 * aJ2 -f aJl + aJp2 
810 E = ref * gain - G 

Z = cl * E + c2 * Epl + c3 * Ep2 + Zpl 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*f<lp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
Zpl = Z 
G = X 

first = 1 
RETURN 

850 ’UVa P I D, equation (4.17), forward path 
IF first = 1 GOTO 860 
K = -359 

LOCATE 3, 22: PRINT ”UVa P I D, equation (4.17)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kp = 1 

Kd = .4 

Kin = .5 

cl = Kp + Kd / T 
c2 = -Kd / T 
dl = T / 2 

860 E = ref * gain - G 

Etotal = Etotal + dl * (E + Epl) 

Z = cl * E + c2 * Epl + Kin * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)-l-al*fd-t-a2*fdpl-(-a3*fdp2-t-a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 

G = X 
RETURN 

900 ’Oak Ridge P I D, equation (4.13), forward path 
IF first = 1 GOTO 910 
K = -100 

LOCATE 3, 19: PRINT ”Oak Ridge P I D, equation (4.13)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kp = 1 

Kdl = .4 

Kd2 = .4 

Kin = .5 * T 

W = .5 

ddO = 1 - W " 3 

ddl = 3 / 2 * (1 - W) ‘ 2 * (1 + W) 
dd2 = 1 / 2 * (1 - W) * 2 
910 E = ref * gain - G 
dE = E - daO 
Etotal = E + Etotal 
daO = daO + ddO * dE 
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dal = dal + ddl * dE 
da2 = da2 + dd2 *dE 

Z = Kp * E + Kdl * dal + 2 * Kd2 * da2+ Kin * Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl+a3*fdp2+a4*fdp3-bUXpl-b2*Xp2-b3*Xp3)/b0 
G = X 
first = 1 
RETURN 

1000 ’Oak Ridge P I D, equation (4.15), forward path 
IF first = 1 GOTO 1010 
K = -400 

LOCATE 3, 19: PRINT "Oak Ridge P I D, equation (4.15)” 

LOCATE 4, 35: PRINT ”K=”; K 

BEEP 

gain = 1 

Kp = 1 

Kpp = 1 

Kdl = .4 

Kd2 = .4 

Kdpl = .4 

Kdp2 = .4 

Kin = .5 * T 

W = .5 

ddO = 1 - W * 3 

ddl = 3 / 2 * (1 - W) * 2 * (1 + W) 
dd2 = 1 / 2 * (1 - W),* 2 
1010 E = ref * gain - G 
dE = E - daO 
Etotal = E + Etotal 
daO = daO + ddO * dE 
dal = dal + ddl * dE 
da2 = da2 + dd2 * dE 
daOp = daO + dal + da2 
dalp = dal +2 * da2 
da2p = da2 

Z = Kpp*da0p+Kp*da0+Kdpl*dalp+2*Kdp2*da2p+Kdl*dal+2*kd2*da2-(-Kin*Etotal 

V = K * Z 

X=(a0*(V+Vpl+Vp2-(-Vp3)-(-al*fd-|-a2*fdpl-|-a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/b0 
G = X 
first = 1 
RETURN 

1100 ’Tustin’s P I D, equation (B-5), forward patli 
IF first = 1 GOTO 1110 
K = -396 

LOCATE 3, 22: PRINT "Tustin’s P I D, equation (B-5)” 

LOCATE 4, 35: PRINT ”K=”; K 
BEEP 
gain = 1 
Kp = 1 


78 



Kd = .4 
Kin = .5 

cl = Kp + 2 * Kd / T + T * Kin / 2 
c2 = T * Kin - 4 * Kd / T 
c3 = T * Kin / 2 + 2 * Kd / T - Kp 
1110 E = ref * gain - G 

Z = cl * E + c2 * Epl + c3 * Ep2 + Zp2 
V = K * Z 

X=(aO*(V+Vpl+Vp2+Vp3)+al*fd+a2*fdpl-l-a3*fdp2+a4*fdp3-bl*Xpl-b2*Xp2-b3*Xp3)/bO 

G = X 

Zp2 = Zpl 

Zpl = Z 

first = 1 

RETURN 



APPENDIX R 


Tustin’s Method of Transformation 

The Tustin’s transformation is a transformation from the 5-domain to the z-domain by 
substituting into the 5-domain equation: 


5 = 


2 0 - 1 ) 

T(z + i) 


, where T is the sampling time 


(B-l) 


The Tustin’s transformation is only an approximation between the 5-domain and z-domain 
which is based on the trapezoidal integration formula. This transformation gives good results as 
long as the sampling rate is high. 

For a dual phase-advance controller given as: 


V _/ l + nA S \ 7 _ 1 -1- 2nA 5 + (nA 5) 2 
f ~\l + A 5 ) ~ i + 2A 5+ (A 5) 2 


then substituting in the Tustin’s transformation of equation B-l, the dual phase-advance has the 
form: 


V k = a ° f k + a l f k-! + ft 2 f k-2 - b l V k-l * b “ V k-2 


(B-3) 


where: 

(T 2 + 4nAT + 4nA 2 ) 
a °“ (T 2 + 4AT + 4A 2 ) 


(2T 2 - 8A 2 ) 

(T 2 + 4 AT + 4A 2 ) 


80 



(2T 2 - 8nA 2 ) _ (T 2 - 4AT + 4A 2 ) 

(T 2 + 4AT + 4A 2 ) 2_ (T 2 + 4AT + 4A 2 ) 

(T 2 - 4nAT + 4nA 2 ) 

(T 2 + 4AT + 4A 2 ) 


For a P I D controller given as: 


V 

c 


K 


[k p 


+ K 


+ S 


(B-4) 


then substituting in the Tustin’s transformation of equation B-4^the P I D has the form: 


V k “ a « f k + a i c k-i + a 2 f k-2 + V k- 


(B-5) 


where: 


2 K i K- T A K , 

a o= ■+ — h ^ a i~ — T~~ 


K- T 2 K, 
a 2 = - K r 


T " P 


The equations B-3 and B-5 arc used as the Tustin’s controller algorithms in the simulation 
program. 

The dual phase-advance controller described by equation B-3 is referred to as Tustin’s 
D P A. The P I D controller described by equation B-5 is referred to as Tustin’s P I D. 
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